
# Estimating the persistence of party cue effects in a panel survey experiment

# PRIMARY ANALYSIS FILE ####

library(brms)
library(tidyverse)
library(tidybayes)
library(cowplot)

set.seed(42)

# Load in the data
df <- read_rds("data_odc.rds")

# Drop people who started but did not complete t0 and t1
df <-
  df %>% 
  filter(finished.t0 == 1 & finished.t1 == 1)

# t0 model -------------------------- ####

# > Specify data for model ####
df_for_t0 <-
  df %>%
  filter(condition == "t0_t1") %>% 
  filter(!is.na(outcome.t0_recode)) # drop missing outcomes

# > Standarize outcome variable ####
control.t0_mean <- mean(df_for_t0[df_for_t0$cue == 0, ]$outcome.t0_recode, na.rm = T)
control.t0_std  <-   sd(df_for_t0[df_for_t0$cue == 0, ]$outcome.t0_recode, na.rm = T)

df_for_t0 <- 
  df_for_t0 %>%
  mutate(outcome.t0_recode_std = (outcome.t0_recode - control.t0_mean) / control.t0_std)

# > Fit model ####

fit_t0 <- 
  brm(data = df_for_t0,
      family = gaussian,
      formula = outcome.t0_recode_std ~ 1 + cue + (1 + cue | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9),
      seed = 42, file = "fit_t0")

summary(fit_t0)


# t1 model -------------------------- ####

# > Specify data for model ####
df_for_t01_vs_t1 <-
  df %>%
  filter(!is.na(outcome.t1_recode)) # drop missing outcomes

# > Standarize outcome variable ####
control.t1_mean <- mean(df_for_t01_vs_t1[df_for_t01_vs_t1$cue == 0 & df_for_t01_vs_t1$t1_only == 0, ]$outcome.t1_recode, na.rm = T)
control.t1_std  <-   sd(df_for_t01_vs_t1[df_for_t01_vs_t1$cue == 0 & df_for_t01_vs_t1$t1_only == 0, ]$outcome.t1_recode, na.rm = T)

df_for_t01_vs_t1 <- 
  df_for_t01_vs_t1 %>%
  mutate(outcome.t1_recode_std = (outcome.t1_recode - control.t1_mean) / control.t1_std)

# > Recode t1_only dummy ####
df_for_t01_vs_t1 <-
  df_for_t01_vs_t1 %>% 
  mutate(t1_only = case_when(t1_only == 0 ~ -0.5,
                             t1_only == 1 ~ 0.5))

# > Fit model ####

fit_t01_vs_t1 <- 
  brm(data = df_for_t01_vs_t1,
      family = gaussian,
      formula = outcome.t1_recode_std ~ 1 + cue*t1_only + (1 + cue*t1_only | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.95),
      seed = 42, file = "fit_t01_vs_t1")

summary(fit_t01_vs_t1)


# Aggregate plots -------------------------- ####

# > t0 data ####
plot_t0 <-
  df %>%
  mutate(cue = factor(cue)) %>% 
  rename(`Party cue` = cue) %>% 
  ggplot(aes(x = condition, y = outcome.t0_recode, color = `Party cue`, shape = `Party cue`)) + 
  theme_bw() +
  geom_point(position = position_jitterdodge(dodge.width = 0.8, jitter.width = 0.5, jitter.height = 0.2), alpha = 0.1) +
  stat_summary(position = position_dodge(0.8), fun.data = "mean_cl_boot") +
  theme(legend.position = c(0.75, 0.95),
        legend.direction = "horizontal",
        legend.box.background = element_rect(),
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text = element_text(color = "black")) +
  scale_y_continuous(labels = 1:7, breaks = 1:7) +
  labs(x = "Condition", y = "Policy opinion (7 = most party-consistent) [95% CI]",
       title = "Data at time 0 (initial survey)")

# > t0-t1 data ####
plot_t01 <-
  df %>%
  mutate(cue = factor(cue)) %>% 
  rename(`Party cue` = cue) %>% 
  ggplot(aes(x = condition, y = outcome.t1_recode, color = `Party cue`, shape = `Party cue`)) + 
  theme_bw() +
  geom_point(position = position_jitterdodge(dodge.width = 0.8, jitter.width = 0.5, jitter.height = 0.2), alpha = 0.1) +
  stat_summary(position = position_dodge(0.8), fun.data = "mean_cl_boot") +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text = element_text(color = "black")) +
  scale_y_continuous(labels = 1:7, breaks = 1:7) +
  labs(x = "Condition", y = "",
       title = "Data at time 1 (follow-up survey 3 days later)")

# > Model estimates plot ####

coef_summary <-
  bind_rows(
  
    # ATE at t0
    posterior_samples(fit_t0, pars = "b_cue") %>% 
      median_hdi(.width = 0.95) %>% 
      rename(value = b_cue) %>% 
      mutate(parameter = "ATE at t0"),
   
    # ATE at t1
    posterior_samples(fit_t01_vs_t1, pars = "b_cue") %>%
      select(b_cue) %>% 
      median_hdi(.width = 0.95) %>% 
      rename(value = b_cue) %>% 
      mutate(parameter = "ATE at t1 (across conditions)"),
    
    # Difference in ATEs at t1 by condition
    posterior_samples(fit_t01_vs_t1, pars = "b_cue:t1_only") %>%
      median_hdi(.width = 0.95) %>% 
      rename(value = `b_cue:t1_only`) %>% 
      mutate(parameter = "ATE at t1: condition t1_only vs. t0_t1")
     
  )

# Plot
plot_coefs <-
  coef_summary %>% 
  ggplot(aes(x = parameter, y = value)) +
  theme_bw() +
  geom_point() +
  geom_errorbar(aes(ymin = .lower, ymax = .upper), width = 0) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_text(aes(label = round(value, 2)), hjust = -0.5) +
  labs(x = "Estimate", y = "Standardized estimate [95% HPDI]",
       title = "Standardized estimates") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text = element_text(color = "black"))

# Join plots
plot_grid(plot_t0, plot_t01, plot_coefs, ncol = 3, labels = "AUTO")

ggsave("fig_2.png", dpi = 300, height = 8, width = 20)


# Issue-level plots -------------------------- ####

# > ATEs at t0 for issues ####
t0_ate_issues <- 
  fit_t0 %>% 
  spread_draws(b_cue, r_item[item, parameter]) %>% 
  ungroup() %>% 
  filter(parameter == "cue") %>% 
  mutate(value = b_cue + r_item) %>% 
  group_by(item) %>% 
  median_hdi(value, .width = 0.95) %>% 
  mutate(parameter = "ATE at t0")

# > ATEs at t1 for issues ####
t1_ate_issues <- 
  fit_t01_vs_t1 %>% 
  spread_draws(b_cue, r_item[item, parameter]) %>% 
  ungroup() %>% 
  filter(parameter == "cue") %>% 
  mutate(value = b_cue + r_item) %>% 
  group_by(item) %>% 
  median_hdi(value, .width = 0.95) %>% 
  mutate(parameter = "ATE at t1 (across conditions)")

# > Difference in ATEs at t1 by condition for issues ####
t1_ate_issues_by_condition <- 
  fit_t01_vs_t1 %>% 
  spread_draws(`b_cue:t1_only`, r_item[item, parameter]) %>% 
  ungroup() %>% 
  filter(parameter == "cue:t1_only") %>% 
  mutate(value = `b_cue:t1_only` + r_item) %>% 
  group_by(item) %>% 
  median_hdi(value, .width = 0.95) %>% 
  mutate(parameter = "ATE at t1: condition t1_only vs. t0_t1")

# Join coefs
coef_issues_summary <- 
  bind_rows(
  t0_ate_issues,
  t1_ate_issues,
  t1_ate_issues_by_condition
  ) %>% 
  # add item indexes from main data
  left_join(.,
            df %>%
              select(item, item_index) %>% 
              distinct(item, .keep_all = TRUE),
            by = "item")

# Plot coefs
coef_issues_summary %>% 
  ggplot(aes(x = item_index, y = value, color = parameter, shape = parameter)) +
  theme_bw() +
  geom_point(position = position_dodge(.9)) +
  geom_errorbar(aes(ymin = .lower, ymax = .upper), width = 0, position = position_dodge(.9)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_text(aes(label = sprintf("%0.2f",round(value, 2))), hjust = -0.5, position = position_dodge(.9), show.legend = F) +
  labs(x = "", y = "Standardized estimate [95% HPDI]",
       title = "Standardized estimates for individual policy issues") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text = element_text(color = "black"),
        legend.position = c(0.5, 0.95),
        legend.box.background = element_rect(),
        legend.title = element_blank(),
        legend.direction = "horizontal")

ggsave("fig_3.png", dpi = 300, height = 8, width = 10)


# Appendix -------------------------- ####

# Model tabulated output and diagnostics ####

# Write function to get table of output
fun_get_table <- function(fit = NULL) {
  
  fit_summary <- summary(fit) # summarise model
  fit_groups  <- unique(fit$ranef$group) # get random fx group names
  
  # Get fixed effects and residual
  tab <-
    bind_rows(
      as_tibble(fit_summary$fixed)     %>% mutate(Term = row.names(fit_summary$fixed),     Group = "fixed"),
      as_tibble(fit_summary$spec_pars) %>% mutate(Term = row.names(fit_summary$spec_pars), Group = "residual")
    )
  
  # Add random fx parameters
  for (i in 1:length(fit_groups)) {
    
    tab <- 
      bind_rows(
        tab,
        as_tibble(fit_summary$random[[i]]) %>% mutate(Term = row.names(fit_summary$random[[i]]), Group = fit_groups[i])
      )
    
  }
  
  tab %>%
    select(Group, Term, everything()) # return output
  
}

# Write function to make traceplots for parameters we care about
fun_get_traceplots <- function(fit = NULL, fit_table = NULL) {
  
  # Get parameter names
  params <- parnames(fit)[1:nrow(fit_table)]
  
  # Get samples
  samples <-
    posterior_samples(fit, pars = params, add_chain = TRUE) %>%
    select(chain, iter, everything())
  
  # Get plots
  map(params, 
      function(.x, samples) {
        
        samples %>%
          rename(value = !!as.name(.x)) %>%
          ggplot(aes(x = iter - 1000, y = value, color = as.factor(chain))) + 
          theme_bw() +
          geom_line(alpha = 0.8) +
          labs(title = paste0(.x), x = "iteration") +
          theme(legend.position = "none",
                plot.title = element_text(hjust = 0.5, face = "bold", size = 8),
                axis.text = element_text(color = "black"))
        
      },
      samples)
  
}

# > t0 model ####

# Table A7
fit_t0_table <- fun_get_table(fit = fit_t0)
saveRDS(fit_t0_table, "table_A7.rds") # write to file

# Figure A2
traceplots_t0 <- fun_get_traceplots(fit = fit_t0, fit_t0_table)
plot_grid(plotlist = traceplots_t0)
ggsave("fig_A2.png", height = 8, width = 8, dpi = 300) # write to file

# > t1 model ####

# Table A8
fit_t1_table <- fun_get_table(fit = fit_t01_vs_t1)
saveRDS(fit_t1_table, "table_A8.rds") # write to file

# Figure A3
traceplots_t1 <- fun_get_traceplots(fit = fit_t01_vs_t1, fit_table = fit_t1_table)
plot_grid(plotlist = traceplots_t1)
ggsave("fig_A3.png", height = 10, width = 12, dpi = 300) # write to file

# Robustness-to-attrition models ####

# > t0 model ####

# >> Adjust for pre-treatment covariates ####

# Mean-center covariates
df_for_t0 <-
  df_for_t0 %>% 
  mutate_at(vars(c(female, college, democrat, party_strength, liberal, ideo_strength, biden2020)),
            .funs = list("c" = ~as.numeric(scale(., center = T, scale = F))))

# Fit model
fit_t0_covariates <- 
  brm(data = df_for_t0,
      family = gaussian,
      formula = outcome.t0_recode_std ~ 1 + cue + 
        female_c + college_c + democrat_c + party_strength_c + liberal_c + ideo_strength_c + biden2020_c +
        (1 + cue | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9),
      seed = 42, file = "fit_t0_covariates")

# Table A9
fit_t0_covariates_table <- fun_get_table(fit = fit_t0_covariates)
saveRDS(fit_t0_covariates_table, "table_A9.rds") # write to file

# Figure A4
traceplots_t0_covs <- fun_get_traceplots(fit = fit_t0_covariates, fit_table = fit_t0_covariates_table)
plot_grid(plotlist = traceplots_t0_covs)
ggsave("fig_A4.png", height = 10, width = 12, dpi = 300) # write to file

# >> Exclude foreign aid issue ####

# Exclude 
df_temp_t0 <-
  df_for_t0 %>% 
  filter(item_index != "Foreign aid")

# Fit model
fit_t0_exclude_fa <- 
  brm(data = df_temp_t0,
      family = gaussian,
      formula = outcome.t0_recode_std ~ 1 + cue + (1 + cue | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9),
      seed = 42, file = "fit_t0_exclude_foreign_aid")

# Table A10
fit_t0_exclude_fa_table <- fun_get_table(fit = fit_t0_exclude_fa)
saveRDS(fit_t0_exclude_fa_table, "table_A10.rds") # write to file

# Figure A5
traceplots_t0_exclude_fa <- fun_get_traceplots(fit = fit_t0_exclude_fa, fit_table = fit_t0_exclude_fa_table)
plot_grid(plotlist = traceplots_t0_exclude_fa)
ggsave("fig_A5.png", height = 8, width = 8, dpi = 300) # write to file

# > t1 model ####

# >> Adjust for pre-treatment covariates ####

# Mean-center covariates
df_for_t01_vs_t1 <-
  df_for_t01_vs_t1 %>% 
  mutate_at(vars(c(female, college, democrat, party_strength, liberal, ideo_strength, biden2020)),
            .funs = list("c" = ~as.numeric(scale(., center = T, scale = F))))

# Fit model
fit_t1_covariates <- 
  brm(data = df_for_t01_vs_t1,
      family = gaussian,
      formula = outcome.t1_recode_std ~ 1 + cue*t1_only + 
        female_c + college_c + democrat_c + party_strength_c + liberal_c + ideo_strength_c + biden2020_c +
        (1 + cue*t1_only | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 42, file = "fit_t1_covariates")

# Table A11
fit_t1_covariates_table <- fun_get_table(fit = fit_t1_covariates)
saveRDS(fit_t1_covariates_table, "table_A11.rds") # write to file

# Figure A6
traceplots_t1_covs <- fun_get_traceplots(fit = fit_t1_covariates, fit_table = fit_t1_covariates_table)
plot_grid(plotlist = traceplots_t1_covs)
ggsave("fig_A6.png", height = 10, width = 12, dpi = 300) # write to file

# >> Exclude foreign aid issue ####

# Exclude 
df_temp_t1 <-
  df_for_t01_vs_t1 %>% 
  filter(item_index != "Foreign aid")

# Fit model
fit_t1_exclude_fa <- 
  brm(data = df_temp_t1,
      family = gaussian,
      formula = outcome.t1_recode_std ~ 1 + cue*t1_only + (1 + cue*t1_only | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 42, file = "fit_t1_exclude_foreign_aid")

# Table A12
fit_t1_exclude_fa_table <- fun_get_table(fit = fit_t1_exclude_fa)
saveRDS(fit_t1_exclude_fa_table, "table_A12.rds") # write to file

# Figure A7
traceplots_t1_exclude_fa <- fun_get_traceplots(fit = fit_t1_exclude_fa, fit_table = fit_t1_exclude_fa_table)
plot_grid(plotlist = traceplots_t1_exclude_fa)
ggsave("fig_A7.png", height = 10, width = 12, dpi = 300) # write to file

# Models excluding Independents ####

# > t0 model ####
fit_t0_excl_inds <- 
  brm(data = df_for_t0 %>% filter(pol_party != 4),
      family = gaussian,
      formula = outcome.t0_recode_std ~ 1 + cue + (1 + cue | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 42, file = "fit_t0_excl_inds")

# Table A13
fit_t0_excl_inds_table <- fun_get_table(fit = fit_t0_excl_inds)
saveRDS(fit_t0_excl_inds_table, "table_A13.rds") # write to file

# Figure A8
traceplots_t0_excl_inds <- fun_get_traceplots(fit = fit_t0_excl_inds, fit_table = fit_t0_excl_inds_table)
plot_grid(plotlist = traceplots_t0_excl_inds)
ggsave("fig_A8.png", height = 8, width = 8, dpi = 300) # write to file

# > t1 model ####
fit_t01_vs_t1_excl_inds <- 
  brm(data = df_for_t01_vs_t1 %>% filter(pol_party != 4),
      family = gaussian,
      formula = outcome.t1_recode_std ~ 1 + cue*t1_only + (1 + cue*t1_only | item) + (1 + cue | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 42, file = "fit_t01_vs_t1_excl_inds")

# Table A14
fit_t1_excl_inds_table <- fun_get_table(fit = fit_t01_vs_t1_excl_inds)
saveRDS(fit_t1_excl_inds_table, "table_A14.rds") # write to file

# Figure A9
traceplots_t1_excl_inds <- fun_get_traceplots(fit = fit_t01_vs_t1_excl_inds, fit_table = fit_t1_excl_inds_table)
plot_grid(plotlist = traceplots_t1_excl_inds)
ggsave("fig_A9.png", height = 10, width = 12, dpi = 300) # write to file

# Models separated by Democrats and Republicans ####

# > t0 models ####

party_names <- as.list(unique(df_for_t0$party_bin))

party_fits_t0 <-
  map(party_names, # maps over Democrat and Republican
      function(.x) {
        
        brm(data = df_for_t0 %>% filter(party_bin == .x),
            family = gaussian,
            formula = outcome.t0_recode_std ~ 1 + cue + (1 + cue | item) + (1 + cue | pid),
            prior = c(prior(normal(0, 0.5), class = Intercept),
                      prior(normal(0, 0.5), class = b),
                      prior(exponential(5), class = sd),
                      prior(exponential(1), class = sigma),
                      prior(lkj(4),         class = cor)),
            iter = 3000, warmup = 1000, chains = 4, cores = 4,
            control = list(adapt_delta = 0.9),
            seed = 42, file = paste0("fit_t0_", .x))
        
      })

names(party_fits_t0) <- party_names # Give names to keep track

# > t1 models ####

party_fits_t1 <-
  map(party_names,
      function(.x) {
        
        brm(data = df_for_t01_vs_t1 %>% filter(party_bin == .x),
            family = gaussian,
            formula = outcome.t1_recode_std ~ 1 + cue*t1_only + (1 + cue*t1_only | item) + (1 + cue | pid),
            prior = c(prior(normal(0, 0.5), class = Intercept),
                      prior(normal(0, 0.5), class = b),
                      prior(exponential(5), class = sd),
                      prior(exponential(1), class = sigma),
                      prior(lkj(4),         class = cor)),
            iter = 3000, warmup = 1000, chains = 4, cores = 4,
            control = list(adapt_delta = 0.95),
            seed = 42, file = paste0("fit_t01_vs_t1_", .x))
        
      })

names(party_fits_t1) <- party_names

# Summarise estimates
party_coef_summary <-
  map(party_names,
  function(.x) { 
    
    bind_rows(
    
    # ATE at t0
    posterior_samples(party_fits_t0[[.x]], pars = "b_cue") %>% 
      median_hdi(.width = 0.95) %>% 
      rename(value = b_cue) %>% 
      mutate(parameter = "ATE at t0"),
    
    # ATE at t1
    posterior_samples(party_fits_t1[[.x]], pars = "b_cue") %>%
      select(b_cue) %>% 
      median_hdi(.width = 0.95) %>% 
      rename(value = b_cue) %>% 
      mutate(parameter = "ATE at t1 (across conditions)"),
    
    # Difference in ATEs at t1 by condition
    posterior_samples(party_fits_t1[[.x]], pars = "b_cue:t1_only") %>%
      median_hdi(.width = 0.95) %>% 
      rename(value = `b_cue:t1_only`) %>% 
      mutate(parameter = "ATE at t1: condition t1_only vs. t0_t1")
    
    ) %>% mutate(party = .x) }
  ) %>% 
  bind_rows()

# Plot
party_coef_summary %>% 
  ggplot(aes(x = parameter, y = value, color = party, shape = party)) +
  theme_bw() +
  geom_point(position = position_dodge(.6)) +
  geom_errorbar(aes(ymin = .lower, ymax = .upper), width = 0, position = position_dodge(.6)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_text(aes(label = round(value, 2)), hjust = -0.5, position = position_dodge(.6), show.legend = F) +
  labs(x = "Estimate", y = "Standardized estimate [95% HPDI]",
       title = "Standardized estimates by party ID") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        axis.text = element_text(color = "black"),
        legend.position = "top",
        legend.title = element_blank()) +
  scale_color_manual(values = c("blue", "red"))

ggsave("fig_A10.png", height = 6, width = 8, dpi = 300) # write to file


# Estimate persistence ratio ----

df_ratio <- 
  df %>% 
  select(pid, item, cue, condition, contains("recode")) %>% 
  pivot_longer(cols = c(contains("recode")),
               names_to = "wave") %>% 
  mutate(wave = str_remove(wave, "outcome."),
         wave = str_remove(wave, "_recode"),
         wave_t1 = case_when(wave == "t0" ~ 0, wave == "t1" ~ 1)) %>% 
  drop_na(value)

# Standarize outcome variable
control.t0_mean <- mean(df_ratio[df_ratio$cue == 0 & df_ratio$wave_t1 == 0, ]$value, na.rm = T)
control.t0_std  <-   sd(df_ratio[df_ratio$cue == 0 & df_ratio$wave_t1 == 0, ]$value, na.rm = T)

df_ratio <- 
  df_ratio %>%
  mutate(value_std = (value - control.t0_mean) / control.t0_std)

# Fit model
fit_ratio <-
  brm(data = df_ratio,
      family = gaussian,
      formula = value_std ~ 1 + cue*wave_t1 + (1 + cue*wave_t1 | item) + (1 + cue*wave_t1 | pid),
      prior = c(prior(normal(0, 0.5), class = Intercept),
                prior(normal(0, 0.5), class = b),
                prior(exponential(5), class = sd),
                prior(exponential(1), class = sigma),
                prior(lkj(4),         class = cor)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 42)

summary(fit_ratio)

# Compute and summarise persistence ratio

p_ratio_samples <- posterior_samples(fit_ratio, pars = c("b_cue", "b_cue:wave_t1"))
saveRDS(p_ratio_samples, "p_ratio_samples.rds")

p_ratio <- 
  readRDS("p_ratio_samples.rds") %>% 
  mutate(cue_at_t1 = b_cue + `b_cue:wave_t1`) %>% 
  mutate(p_ratio = cue_at_t1 / b_cue) %>% 
  pivot_longer(cols = everything(),
               names_to = "parameter") %>% 
  group_by(parameter) %>% 
  median_hdi(.width = 0.95)

# Table A15
fit_ratio_table <- fun_get_table(fit = fit_ratio)
saveRDS(fit_ratio_table, "table_A15.rds") # write to file
